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Abstract 


This  is  a  comparison  study  of  the  abilities  of  the  eigenvalue  method 
as  a  numerical  method  in  solving  the  transient  heat  conduction  equation. 
The  eigenvalue  method  was  compared  to  five  other  numerical  methods: 
Runge-Kutta,  Gears,  extrapolation,  fully  implicit,  and  Crank-Nicolson. 
These  methods  were  used  to  solve  three  physical  problems.  The  first  is 
a  two  dimensional  slab  which  takes  advantage  of  the  symmetry  of  the 
problem.  The  second  is  the  same  slab  problem  without  taking  advantage  of 
the  symmetry.  And  the  third  is  a  cylindrical  problem  taking  full  advan¬ 
tage  of  symmetry. 

The  scope  of  the  study  is  to  see  which  methods  take  less  computer 
time  while  maintaining  sufficient  accuracy.  The  time  it  takes  the  com¬ 
puter  to  totally  execute  the  program  was  used  as  the  time  comparison 
basis.  The  accuracy  is  a  comparison  of  the  exact  solution  to  the  numeri¬ 
cal  solution.  A  root  mean  square  average  of  all  the  grid  points  per  time 
step  is  used. 

The  results  of  the  study  were  surprising.  The  accuracy  of  the 
eigenvalue  method  is  not  any  better  than  that  of  the  Crank-Nicolson 
method.  The  computer  times  show  that  the  eigenvalue  is  not  the  fastest 
for  short  transient  times.  A  long  transient  problem  with  nonlinear  terms 
was  not  used  in  this  study. 


A  COMPARISON  STUDY  OF  THE  EIGENVALUE  METHOD  FOR 


THE  SOLUTION  OF  THE  TRANSIENT  HEAT  CONDUCTION  EQUATION 

I .  Introduction 

In  the  field  of  nuclear  engineering,  the  study  of  reactor  thermal- 
hydraulics  is  coupled  to  the  study  of  reactor  neutronics.  This  is  because 
reactor  power  is  directly  proportional  to  neutron  population.  Active 
core  volumes  are  getting  greater  than  300,000  liters  while  the  neutronics 
people  are  calculating  the  multiplication  factor  to  four  decimal  place 
accuracy  (8:634). 

Solving  these  problems  takes  massive  computer  codes,  such  as  COBRA, 
to  solve  just  the  thermalhydraulic  problem.  Solution  of  the  transient 
heat  conduction  equation  is  only  a  small  part  of  the  massive  problem. 

By  solving  parts  of  the  big  problem  faster,  the  overall  solution  could  be 
done  faster  and  more  accurately. 

The  scope  of  the  study  is  limited.  Density,  specific  heat,  and 
thermal  conductivity  have  all  been  assumed  constant.  These  were  then 
lumped  into  the  single  constant,  thermal  diffusivity.  This  makes  the 
solution  limited  but  it  allows  the  use  of  dimensionless  parameters.  None 
of  the  problems  considered  had  a  heat  source.  The  temperature  difference 
used  was  200°  K. 

The  problems  were  investigated  using  Cartesian  and  cylindrical 
coordinates.  Three  problems  were  looked  at.  Two  problems  were  actually 
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one  problem  solved  two  different  ways.  One  way  took  advantage  of  symmetry 


and  the  other  way  did  not.  The  problem  was  a  two  dimensional  slab  with 
Dirichlet  boundary  conditions.  The  third  problem  was  a  cylindrical  problem 
with  Dirichlet  boundary  conditions.  Full  advantage  was  taken  of  the 
symmetry  in  the  cylindrical  problem. 

Six  different  numerical  schemes  were  used  to  solve  the  problems: 
eigenvalue,  Runge-Kutta,  extrapolation.  Gears,  implicit,  and  Crank-Nicolson 
methods.  The  computer  used  was  AFIT's  SSC  VAX  11/780.  Each  computer  run 
was  timed  for  the  amount  of  CPU  time  and  wall  clock  time  necessary  to 
execute  the  command. 


II .  Theory 

Solution  of  the  transient  heat  conduction  equation  is  important  in 
many  engineering  problems.  It  is  of  even  greater  importance  for  reactor 
design  since  the  amount  of  heat  a  reactor  produces  is  directly  coupled 
to  the  amount  of  heat  removed.  This  study  looked  at  different  numerical 
methods  for  solving  the  transient  heat  conduction  equation: 

=  aV2T  (l) 


Here  T  is  temperature,  t  is  time,  a  is  thermal  diffusivity,  and  V2  is  the 
Laplacian  operator. 

The  transient  heat  conduction  equation  was  solved  by  three  different 
methods:  implicit,  Crank-Nicolson ,  and  the  method  of  lines.  The  theory 

of  each  method  will  be  described  briefly. 

Implicit  Method 

The  implicit  method  (IM)  is  in  common  use  among  engineers.  This  is 
because  of  its  unconditional  stability.  The  method  generates  systems  of 
simultaneous  equations  due  to  the  differencing  equations  used.  The 
differencing  equations  discretize  the  spatial  variables  with  respect  to 
the  next  time  step  using  a  central  difference  (equations  (2)  and  (3)). 

A  standard  forward  difference  is  used  on  the  time  derivative  (equation 
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where  k  is  the  time  step,  i  and  3  are  spatial  nodal  points. 

Plugging  the  above  difference  equations  into  equation  (1)  and 

solving  for  T  : 
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where  (A  x  )  2  =  (A  y )2  =  h2  and  C  =  a  ^  ; t *  .  The  above  equation 

produces  a  set  of  simultaneous  equations  with  errors  of 

0(  h  )  +  0(  h  ) 2 . 
t  x ,  y 


C rank -Nico Ison 

A  method  with  improved  accuracy  and  unconditional  stability  is  the 
Crank-Nicolson  (CN)  method.  The  CN  method  uses  an  averaging  scheme  of 
central  time  differences;  see  equations  (6)  and  (7). 
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The  time  derivative  is  a  central  difference  at  the  k  +  time 
node.  Taking  equations  (4),  (6),  and  (7)  and  substituting  them  into 
equation  (1)  and  putting  all  unknown  temperatures  (k+1)  on  the  left  side, 
results  in: 
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.  Again  we  develop  a  set  of  simultaneous  equations. 


The  important  point  is  the  error  has  improved  to  0(h  )*  +  0(h  )i 
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Method  of  Lines 

A  numerical  method,  uncommon  for  solving  partial  differential  equa¬ 
tions,  is  the  methods  of  lines  (MOL).  This  method  changes  a  partial 
differential  equation  into  an  ordinary  differential  equation  (ODE)  of 
the  form: 

=  f(t,y)  (9) 

Numerous  books  are  available  for  methods  of  solving  equation  (9) 

(  1 ; 10 ; 2 0  ) .  A  formal  proof  of  turning  equation  (1)  into  the  form  of 
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equation  (9)  is  in  either  Shih  and  Skladany  or  Landry  (25:411—412; 19: 11-15) 
In  practice,  the  spatial  variables  are  discretized  leaving  the  time 
derivative  alone.  Here,  the  spatial  variables  will  be  discretized  with 
a  central  difference. 


3x2  (Ax)2  (10) 
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After  differencing  the  spatial  variables  and  substituting  into 
equation  (1),  the  resulting  equation  is 
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This  is  an  ODE.  The  three  1MSL  (15)  routines  used  were  DREBS,  DVERK,  and 
DGEARS.  They  use  extrapolation,  Runge-Kutta,  and  predictor-corrector 
methods,  respectively.  All  three  routines  were  designed  to  solve  ODE’s 
of  the  form  of  equation  (9)  with  an  initial  condition  supplied. 

DREBS .  Extrapolation  is  the  method  used  in  this  IMSL  routine. 

The  actual  method  is  a  modified  form  of  a  subroutine,  DESUB,  developed 
by  Bulirsch  and  Stoer  (2:10-13).  It  uses  rational  function  for  the  extrap¬ 
olation.  The  exact  error  figure  is  not  given  in  the  documentation. 

Error  bounds  and  stability  are  still  unanswered  questions  for  extrapolation 


methods  (10:101). 


No  documentation  explaining  the  theory  of  the  method  was  available 


in  English.  Therefore,  only  the  basic  algorithm  used  by  Bulirsch  and 
Stoer  will  be  given  (2:3). 
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DVERK .  The  RK  method  is  one  of  the  most  common  methods  for  the 
solution  of  ordinary  differential  equations.  It  is  a  one  step  method, 
only  requiring  the  information  at  time  step  k  to  get  the  information  at 
time  step  k+1.  Therefore  the  method  is  self  starting.  The  RK  method 
does  not  require  the  evaluation  of  any  derivatives  of  f(t,y). 

An  example  of  a  second  order  RK  method  is  the  improved  Euler  method. 
With  the  help  of  Figure  1,  a  geometric  interpretation  on  how  the  method 


works  can  be  shown.  What  is  done  is  an  averaging  of  slopes.  The  point 

(t^  +  h,  +  hy^)  is  computed  using  Eulers  method  giving  the  slope  of 

L  .  Now  the  slope  of  the  curve  is  computed  at  (t  +  h,  y  +  hy  ),  this 
1  K  K  k 

gives  the  slope  of  L^.  The  slopes  of  L  and  are  averaged  giving  the 
slope  of  L  .  A  line  is  drawn  parallel  to  L,  L.  Where  L  intersects  the 


ordinate  axis  at  t  =  t,  ,  =  t,  +  h  is  the  next  point  (t,y). 

k+1  k 

The  slope  of  L  and  L  is  given  by 


\  =  2  f(V  V  +  f(tk  +  h*  yk  +  hV 
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The  equation  of  line  L  is  given  by: 


y  *  yk  +  (t“V  \ 


therefore , 


yk+l  =  yk  *  h  \ 


Equations  (18),  (19)  and  (21)  define  the  improved  Euler  method,  some¬ 
times  called  Heun's  method  (22:317-320). 

This  method  requires  only  two  evaluations  of  the  function  f(t,y). 


once  at  (t  ,y  )  and  once  at  (t  ♦  h,  y  +  hy  ) .  For  comparison,  a 

K  K  K  K  K 

Taylor  series  of  the  same  order  would  require  three  function  evaluations 

f(t,y),  f  and  f  where  f  and  f  are  the  partial  derivatives  of  f(t,y) 
x  y  x  y 


with  respect  to  the  subscript. 


The  most  common  RK  method  is  a  fourth  order  method  (18:835).  What 
makes  this  method  so  appealing  is  the  fact  that  no  derivatives  need  to 
be  calculated,  the  function  itself  only  needs  to  be  evaluated,  and  the 
method  is  self  starting.  The  method  only  has  partial  stability  though. 
Gears  method,  which  is  a  predictor-corrector  method,  does  have  complete 
stability  and  an  ability  to  handle  stiff  systems  of  differential  equations. 

DGEARS ■  Predictor-corrector  methods  are  another  common  numerical 
technique  for  the  solution  of  ODE'S.  The  major  draw  back  with  these 
methods  is  that  they  are  not  self  starting.  A  method  like  RK  is  needed 
to  start  predictor-corrector  methods.  Once  enough  points  are  known  the 
predictor-corrector  technique  is  employed  for  the  rest  of  the  solution. 

The  reason  predictor-corrector  methods  are  not  self  starting  is  that  they 
require  values  at  prior  nodal  points. 

For  a  brief  explanation  of  the  predictor-corrector  formulas,  a 
second  order  formula  will  be  used.  This  should  show  how  the  method 
works  and  also  give  a  comparison  with  the  RK  method  just  described. 

First  a  predicted  value  of  y^  ^  is  needed.  This  is  done  by  finding 
the  slope  at  (t  ,y  )  and  a  line  L  is  drawn  through  this  point,  as  in 
Figure  2.  Next  draw  a  line  L  parallel  to  L  but  through  point  (t 

1  K  -  i  , 

^k  1*'  where  this  line  L,  crosses  the  value  of  t  ^  gives  the  first  guess 
( 0 ) 

of  (t  ,,y  ,).  The  superscript  here  (0)  indicates  that  this  is  the 

k-t-1  k  + 1 

first  guess  at  y^_  --a  predicted  value.  The  equation  used  to  get 
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The  ^  used  in  the  above  equation  is  why  this  method  is  not  self 
starting . 

(  0  ) 

Now  it  is  time  to  correct  the  first  value  of  y^  .  A  slope  can  be 

calculated  at  (t,  , , y,  ,)  which  gives  us  line  L  in  Figure  3.  Lines  L 

k+1  k+1  2  1 

and  are  then  averaged,  getting  line  L.  A  line  L  is  then  drawn  parallel 

to  L  through  the  point  (t,  , y,  ).  The  intersection  at  t  =  t  ,  will  give 

k  k  k  +  1 

the  corrected  value  of  y,  ,  or  y,^  ] .  With  this  new  y^ ^  a  better  esti- 

k+1  k+1  k+1 

mate  of  the  slope  of  is  gotten.  Using  this  better  approximation  of 
the  slope  of  L  a  new  estimate  of  y  is  gotten  or  y^2'.  In  general  it 

Z  K+1  K+1 

is  given  by 
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The  iterations  are  stopped  when 
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where  e  is  the  convergence  criteria  of  the  user.  It  can  be  shown  that 
the  corrector  formula  does  converge  (22:332-333).  If  more  than  two 
iterations  of  equation  (23)  are  required  to  meet  the  convergence  criteria. 


then  the  step  size  is  too  big  and  should  be  reduced. 
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Equations  (22)  and  (23)  define  a  second  order  predictor-corrector 
method.  There  are  obvious  similarities  between  the  RK  and  predictor- 
corrector  methods. 

Eigenvalue  Method .  The  eigenvalue  method  (EM)  should  fall  under 
the  MOL  for  it  also  turns  a  partial  differential  equation  into  an  ODE. 

A  formal  proof,  as  stated  earlier,  of  why  the  EM  works  may  be  found  in 
either  Shih  and  Skladany  or  Landry. 

Taking  mainly  from  Landry:  once  the  differential  equation  is  in 
the  matrix  form, 

{ e }  =  [D]  {9}  (25) 

Where  {}  is  used  for  an  n  component  vector  and  [  ]  are  used  for  an 

n  x  n  matrix.  This  is  the  matrix  form  of  equation  (9).  Where  the 
solution  has  the  form  of 

{ 0 ( t )  }  =  l{c.}  exp({Xi}t)  (26) 

i 

{ c ^ }  are  undetermined  constants  and  fX.}  are  the  eigenvalues  of  the 
system  of  equations.  The  constants  fc  1  turn  out  to  be  another  constant 

1  l  ‘ 

{a^}  multiplied  by  the  eigenvector  { v ^ .  Therefore,  the  final  solution 
in  matrix  form  is 


{t,}  =  {Tf}  -  l  {a.}  {v?}  exp  ({X.}  t)  (27) 

i  .  3 

where  are  the  boundary  conditons,  v^  is  the  jth  component  of  the  ith 
eigenvector  associated  with  the  ith  eigenvalue.  Expanded  the  equation 


looks  like 


III.  Applications 


More  emphasis  was  placed  on  solving  the  transient  heat  conduction 
equation  by  the  MOL  than  by  the  IM  or  CN  methods.  Landry  in  his  study 
and  also  Shih  and  Skladany  in  theirs  did  not  put  much  emphasis  on  the 
MOL  technique. 

The  basic  problem  looked  at  is  a  piece  of  material  (slab  or  cylinder) 
at  some  initial  temperature,  T0 ,  suddenly  having  its  outside  boundaries 
raised  instantaneously  to  a  temperature  of  200  degrees  Kelvin  above  the 
initial  temperature,  T0 .  Three  problems  were  looked  at:  first,  a  slab 
was  analyzed  taking  full  advantage  of  symmetry;  secondly,  the  same  slab 
was  analyzed  without  taking  advantage  of  symmetry;  and  thirdly,  a 
cylinder  was  analyzed  taking  full  advantage  of  symmetry. 

Error  Analysis 

To  help  evaluate  the  method's  accuracy,  the  exact  solution  was  also 
found  for  each  particular  problem.  All  the  grid  points,  for  a  given 


time  s 

tep,  sho 

uld 

be  considered 

to  get 

an  overall  error 

(27: 

:  1023  ) 

The 

exact 

solution 

is 

used  along  with  the 

numerical  solution 

at 

each 

grid 

point 

to  get  a 

root  mean  square 

(  RMS ) 

average  value  for 

the 

time 

step . 

where  u^  is  the  value  of  the  temperature  at  the  ith  nodal  point,  u  is 
the  exact  solution  at  that  point,  and  N  is  the  total  number  of  grid 


points . 


In  the  exact  solution,  only  20  terms  were  used.  This  is  because 


about  15  terms  are  needed  to  get  it  to  converge  to  10  accuracy 

(19:27).  This  accuracy  is  more  than  sufficient  for  engineering  applica- 
-3  -4 

tions  where  10  to  10  accuracy  is  usually  the  norm. 


Symmetric  Slab  Problem 

Figure  4  shows  the  grid  system  used.  It  has  already  been  noted  in 
Landry  that  the  center  of  the  slab  should  have  a  node  also  (19:18). 
Obviously,  it  is  a  crude  assumption  that  the  center  of  the  slab  is  at 
the  same  temperature  as  node  20.  It  is  hard  to  see  why  a  node  21  was 
not  used.  The  original  authors  were  taking  advantage  of  the  symmetry  in 
the  problem. 

The  exact  solution  of  this  problem  is 

r  r  16  /n-nx  \  /miTy\ 

0(x-y't)  =  l  l  ^  sin  (t)  sin  hr) exp 

n  m 

(  30) 

The  IM  method  used  equation  (5)  directly,  and  where  i  is  the  x 

spatial  nodes  and  j  is  the  y  spatial  nodes.  The  CN  method  did  not  use 

k  1 

equation  (8)  directly.  The  CN  method  solved  equation  (8)  for  T  first 


rk+ 1 

1.3 


(1-2C)  Tk  C 

(  1  +  2C  )  i  ,  3  2  (  1  +  2C  ) 


r 


T 


k  +  1 

i+  1 ,  j  + 


k  + 1 

ri-  1  , 


+  T 


k  +  1 
i ,  3  +  1 


+  T 


k+  1 
1.3-1 


k  k  k  k 

+  T  +  T  +  T  +  T 

i+l,j  1-1.3  1,3+1  1,3-1 


(  31  ) 


This  equation  was  used  in  programing  the  CN  method. 
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The  MOL  uses  equation  (12)  directly  for  the  recursion  formula.  The 
difference  between  the  IMSL  routines  and  the  EM  is  in  the  use  of  equation 
(12).  The  IMSL  routines  all  need  to  have  the  first  node  on  the  boundary 
as  the  starting  node.  The  IMSL  routines  all  require  an  external  subroutine 
be  provided.  This  subroutine  evaluates  equation  (12)  at  each  nodal  point 
(see  Appendix).  The  EM  only  used  equation  (12)  to  generate  the  matrix  [D] 
of  equation  ( 25  ) . 

Results .  Two  things  are  being  looked  at  in  the  study:  first  is  the 
time  the  program  takes  to  run  and  second,  the  amount  of  error  involved 
with  the  results. 

To  get  the  time  required  by  the  computer  to  execute  the  program  the 
UNIX  time  command  was  used  (23:473-474).  The  actual  command  used  on 
AFIT's  SSC  computer  was  /bin/time  instead  of  simply  time.  This  returned 
real,  user,  and  system  times  all  in  seconds.  Real  time  is  the  wall  clock 
time,  user  and  system  time  are  CPU  time.  The  results  of  the  timed  computer 
runs  can  be  seen  in  Table  I. 

The  symmetric  problem  has  shown  some  different  results  than  those  of 
Landry.  The  wall  clock  time  is  insignificant  because  the  computer  is  a 
time  sharing  machine.  The  CPU  time  should  not  vary  significantly,  for 
this  is  the  amount  of  time  the  computer  spent  working  on  the  problem. 

Landry  shows  a  CPU  time  for  the  EM,  IM,  and  CN  methods  as  0.8,  4.2,  and 
6.4  respectively  (19:29). 

When  the  computer  runs  were  timed,  the  whole  program  was  timed, 
that  is,  the  setting  up  of  variables,  reading  in  data  from  files,  and 
writing  the  results  to  files.  This  could  explain  why  the  results  of  the 
time  runs  are  so  much  greater  than  those  of  Landry. 
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TABLE  I 

Time  Comparison  of  Symmetric  Slab  Problem 


Method 

CPU  (sec) 

Wall  Clock 

Eigenvalue 

7.1 

19.1 

Runge-Kutta 

0.9 

1.9 

Gears 

1.8 

2.9 

Extrapolation 

11.9 

37 . 7 

Implicit 

50.6 

166.9 

C rank -Ni col son 

50.5 

157.5 

The  errors,  as  compared  in  Table  II,  can  be  seen  to  be  as  good  as 
those  obtained  by  the  EM.  Even  the  ODE  solvers  have  errors  comparable  to 
that  of  the  EM.  These  errors  are  shown  with  respect  to  the  ten  time  steps 
TBAR.  The  closer  the  method  is  to  the  exact  solution,  the  smaller  the 
value  of  the  error. 

Landry  has  shown  that  a  comparison  of  the  accuracies  of  the  EM  with 

the  exact  solution  has  errors  of  ±0.3%  after  the  second  time  step  (19:23). 

Before  this  time  step  errors,  as  compared  with  the  exact  solution,  are 

considerably  worse.  Shih  and  Skladany  noticed  this  problem  and  came  up 

with  an  explanation  for  the  cause  (25:417). 

Since  initially  the  heating  rate  9T/3t  and  the  change  of  the 
heating  rate  32T/3t2  are  large,  and  since  3T/3t  is  related  to 
the  space  derivatives  by 


it  follows  that  initially  the  value  of  V4T  must  be  large  because 


Ivl 
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Now  we  already  know  that  the  fourth  derivative  V 4T  is  embedded 
in  the  leading  truncation  term  in  the  second-order  accurate 
finite-difference  approximation. 


This  states  that  the  truncation  error  of  the  initial  heating  rate  will 
be  large  in  any  scheme  using  a  central  differencing  approximation  of 


the  spatial  variables. 


Full  Grid  Slab  Problem 

Many  engineering  problems  involve  many  more  nodal  points  than  just 
20.  Typically,  more  than  a  hundred  nodal  points  are  used  in  most  applica¬ 
tions.  To  get  a  better  idea  on  how  these  methods  compare  in  a  big  grid, 
the  original  slab  problem  was  rerun  without  taking  advantage  of  symmetry. 
The  number  of  nodal  points  went  from  20  to  121.  All  the  methods  were 
rerun  for  the  full  grid  case  except  for  the  IM  and  CN  methods.  Problems 
arose  in  solving  the  set  of  simultaneous  equations  and  the  results  were 


not  correct. 


All  of  the  equations  used  for  the  full  grid  case  were  the  same  as  in 


the  symmetric  case.  The  same  IMSL  routine  used  to  solve  the  set  of 
simultaneous  equations  in  the  symmetric  case  was  also  used  in  the  full 
grid  case. 

Results .  Table  III  shows  the  timed  run  results  for  all  the  MOL 
programs.  Even  though  the  IM  and  CN  methods  did  not  run  properly,  from 
a  comparison  of  Landry's  results  of  the  symmetric  and  full  grid  times  and 
the  results  that  this  study  got  for  times,  it  is  felt  that  the  RK  and 
Gears  methods  are  faster  than  the  IM  and  CN  methods. 


TABLE  III 


Time  Comparison  for  Full  Grid  Problem 


Method 

CPU  (sec) 

Wall  Clock  (sec) 

Eigenvalue 

462.2 

659.9 

Runge-Kutta 

6.9 

25.0 

Gears 

8.4 

24.0 

Extrapolation 

71.0 

229.9 

The  error  comparison  was  done  for  the  full  grid  problem  also.  The 
results  obtained  do  not  agree  with  those  of  Landry.  As  can  be  seen  in 
Table  IV,  the  error  for  the  EM  is  comparable  to  all  the  other  methods 
(the  values  for  the  IM  and  CN  methods  came  from  Landry's  thesis). 


Cylindrical  Problem 

The  most  common  nuclear  engineering  geometry  is  a  cylinder.  There¬ 
fore,  most  nuclear  engineering  problems  are  done  in  cylindrical  coordinates 
The  third  and  final  problem  looked  at  was  a  solid  cylinder.  The  differen¬ 
cing  equations  were  all  done  in  cylindrical  coordinates.  The  first 
derivative  of  temperature  with  respect  to  the  radius  was  set  equal  to 
zero  at  the  center  of  the  cylinder.  The  rate  of  change  of  temperature 
with  respect  to  the  angle  0,  was  set  equal  to  zero.  The  Laplacian  of 
this  problem  then  is 


rw 


"■'‘7V 


and  when  looking  at  points  on  the  cylinder  centerline  it  becomes 


V2T 


92T  (  3^T 

3r2  +  9z* 


(33) 


The  two  in  equation  (33)  comes  from  an  expansion  done  on  the  middle 
term  in  equation  (32).  The  value  of  r  at  the  centerline  and  9T/9r  at  the 
centerline  are  both  zero  giving  an  indeterminant  value  there.  From  an 
expansion  of  this  term  it  can  be  shown  that 


lim  1  9T  _  9J_T 
r  -*•  0  r  9r  9r^ 


Therefore  there  is  a  two  in  equation  (33). 

A  central  difference  was  used  on  the  center  term  of  equation  (32) 


1  II 

r  3 r 


_k  k  k  k 

T  -  T  T  -  T 

1  i»l.j  =  i  +  1.3  i-l.j 

i(Ar)  2(Ar)  “  2i(Ar)i 


(  34  ) 


The  recursion  formula  used  in  the  cylindrical  problem  using  the  XM 
method  was 


T1  .  =  (  1  +  4C)  Tk  +  1  -  C 

k.3  1.3 


(l  *  -M  Tk  +  1  ♦  (l  -  —  )  Tk  +  1 

\  2 i  /  i  +  1. 3  \  2 i  /  i-1, j 


k+1  k+1 

¥  T  +  T 

i.  3  +  1  i. 3-1 


(  35) 


Then  for  the  centerline 
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The  formula  used  in  the  CN  method  was 


( 1+2C )  Tk+1  -  \  C 

1.1  2 


(l*  if) 


mk+l 

T  .  .  + 


i  +  1.  j 


(-  if) 


+  T 


k+1 
i.  j-!_ 


k 

1 

/  1  \ 

=  (  1-2C )  T. 

+  -  C 

(1+  —  ) 

i .  1 

2 

\  4 1  / 

k  k 

+  T  +  T . 
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Then  when  evaluating  the  problem  on  the  centerline 


,,  _k+l  1  „  ,mktl  k+1  k+1  k  + 

(  1  +  3C  )  T  --C2T  +2T  +T  +T 

i.l  2  L  1+1>1  i-l.j  1,1+1  i. 


1  —  3C  )  T 


k  1  ,  k  k  k 

.  .  +  -  C  2T .  .  +  2T.  ,  .  +  T. 

i.l  2  1+1,1  i-l.l  i.l+l 
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The  basic  formula  to  solve  the  problem  by  the  MOL  is 


~  I 

dt  i  ,  j 


T .  ,  .  -  2T.  .  +  T .  , 

1  +  1.1 _ i.l  i-l.l 


i  +  l.l 


(Ar)' 


+  2  i  ( A  r  ) 


T.  ,  -  2T  +  T 

1.1  +  1 _ i.l  b 


(  Az  )  ' 


and  when  evaluating  on  the  centerline 


Full  advantage  was  taken  of  symmetry.  Figure  5  shows  the  grid 
scheme  used  in  this  study.  When  i=5  ,  then  equations  (36),  (39),  and 

(40)  were  used. 

The  exact  solution  of  this  problem  used  was  (21:165) 


0  (  r ,  z  ,  t : 
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where 
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and  p  are  the  roots  of  the  Bessel  function  of  the  first  kind  of  zero 
n 

order . 

This  exact  solution  was  found  to  have  some  problems.  The  values  of 
the  summation  should  be  between  zero  and  one.  Some  of  the  values  were 


greater  than  one  causing  erroneous  values  in  the  solution  (i.e.,  negative 
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TABLE  V 


Time  Comparison  for  Cylindrical  Problem 


Wall  Clock  ( sec  ) 


Runge-Kutta 


Extrapolation 
Implicit  Method 


C rank -Nico Ison 


temperatures).  All  negative  temperatures  were  set  to  T0  assuming  no 
temperature  change  at  that  nodal  point. 

Result .  The  computer  time  needed  to  run  this  problem  is  found  in 
Table  V.  As  before,  the  RK  method  was  the  fastest. 

The  error  comparison  can  be  found  in  Table  VI.  The  results  shown 
here  are  consistent  with  the  results  for  the  previous  two  problems. 

That  is,  the  EM  does  not  show  to  have  significantly  better  accuracy  than 
the  other  methods  (19:31). 

The  explanation  for  this  is  in  the  differencing  scheme.  The  partial 
differentia',  equation  was  turned  into  an  ODE  by  the  use  of  a  central 
difference  on  the  spatial  variables.  That  central  difference  has  trunca¬ 
tion  error  of  0(h)2.  The  IMSL  subroutines  demonstrate  this  fact  in  the 
accuracies  being  almost  identical  with  the  EM  for  all  three  problems 
examined . 
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IV.  Conclusions  and  Recommendations 


Conclusions 

The  EM  is  not  as  accurate  as  Landry  had  said  (19:44).  The  results 
obtained  have  indicated  that  the  EM  is  as  accurate  as  the  CN  or  IM 
methods.  This  is  believed  to  be  connected  to  the  Taylor  series  differen¬ 
cing  used  on  the  spatial  variables.  The  results  have  shown  that  the  ODE 
solvers  used  from  IMSL  are  as  accurate  as  the  EM  even  though  the  IMSL 
routines  use  other  methods  with  truncation  errors  to  solve  the  ODE. 

The  EM  is  probably  the  easiest  method  to  program.  Once  the  [D] 
matrix  of  equation  (25)  is  set  up,  it  can  be  loaded  into  the  program  via 
an  input  file.  Once  the  eigenvalues  and  eigenvectors  are  found,  the 
number  of  time  steps  that  is  taken  can  be  changed  to  whatever  is  needed 
or  even  made  an  input  variable. 

As  far  as  the  fastest  routines,  the  EM  is  not  as  fast  as  the  RK 
method.  The  RK  method  is  the  fastest  but  not  as  easy  to  program.  The 
problems  run  here  have  the  evaluation  of  the  derivative  being  very 
inexpensive  in  computer  time.  This  tends  to  make  the  RK  method  better 
(14:617). 

Overall  the  MOL  is  the  better  method  over  CN  and  IM.  The  programs 
are  easier  to  write  because  a  canned  subroutine  is  available  to  do  all 
the  work.  No  less  accuracy  was  encountered  when  using  the  IMSL  ODE 
solvers . 

Recommendat ions 

One  thing  that  has  not  been  done  yet  is  to  use  the  EM  on  a  nonlinear 


problem.  An  example  would  be  to  have  heat  loss  by  radiation.  The 
example  in  Shih  and  Skladany  required  the  eigenvalues  and  eigenvectors 
to  be  found  each  and  every  time  step.  Finding  eigenvalues  and  eigen¬ 
vectors  is  another  area  being  looked  into  heavily  (19:46).  For  all  the 
problems  presented  here  the  eigenvalues  and  eigenvectors  were  found 
only  once  for  all  ten  time  steps. 

Another  thing  would  be  to  time  the  EM  for  long  transient  times 
(i.e.,  50  or  more  time  steps),  then  compare  the  amount  of  computer  time 
needed  with  that  of  the  other  MOL  techniques. 

If  the  problem  does  not  have  any  nonlinear  terms,  then  the  eigen¬ 
values  and  eigenvectors  only  need  to  be  found  once.  Where  the  MOL 
requires  the  evaluation  of  the  derivatives  more  than  once  per  time  step, 
some  interesting  results  could  be  derived  from  this  test. 

Finally,  a  problem  more  like  those  found  in  nuclear  engineering, 
consider  a  cylinder  with  a  heat  source  in  the  center  surrounded  by 
cladding.  The  first  problem  here  is  to  find  the  exact  solution.  This 
is  a  coupled  partial  differential  equation,  which  solving  for  the  exact 
solution  cannot  be  done  by  standard  methods. 


Appendix:  Subroutines  Used  in  the  Study 

This  appendix  contains  the  subroutines  used  in  the  IMSL  routines 
for  this  study.  This  is  provided  for  the  reader  to  help  them  better 
evaluate  the  routines  used. 

Any  questions  concerning  the  IMSL  routines  should  be  directed  to 


the  company  given  in  the  bibliography. 


SUBROUTINE  FCN1  (N, X, Y, YPRIME) 

INTEGER  N 

PEAL  Y  ( N  >  , YPR I ME ( N ) , X 
YPRIME  ( 1 )  -  0.0 

YPR  I  ME  ( 2 )  ■  ((2*Y(1))  -  (4*Y(2>>  «■  (2*Y(3>>) 

YPRIME  !  3 )  -  (Y<1)  ♦  Y(2>  +  Y(4>  ♦  YO)  -  (4#Y(3>>> 

YPRIME 14)  -  <  Y ( 1 )  +  Y  <  3)  +  Y<3)  *  Y<9>  -  (4*Y(4)>) 

YPRIME  ( 3)  -  (Y'l)  ♦  Y 14)  *  Y(6)  +  Ytl0>  -  <4*Y(3)>) 

YPRIME  ( 6)  -  (Y<1)  ♦  Y<3>  «■  Y(7>  YII1)  -  <4#Y(6>)> 

YPRIME ( 7)  *  (Y<1)  ♦  (2*Y(6>>  ♦  YI12I  -  (4*Y(7>)> 

YPRIME  <G)  -  ('.2*Y<3)>  ♦  (2»Y(9)>  -  <4*Y18>>> 

YPRIME  ( 9 )  -  ( Y  <  4 )  ♦  YO)  +  Y(10>  +  Y113I  -  (4*Y(9>>> 
YPRIME  (  10)  =*  ( Y !  3 )  ♦  Y<9)  «■  YIU)  +  Y  (  14)  -  (41YI10D) 
YORIMEUl)  «  <  Y  ( 6 )  ♦  Y(10>  +  Y  ( 12)  ♦  YUS)  -(<t*Y(lllll 
YPRIME  (12)  »  (Y  (7)  ♦  I21YU1))  «•  Y  (  16)  -  (4*YI12)I) 

YPRIME  (13)  ■  <  ( 2#Y  <  9 )  )  ♦  (2*YIM)I  -  (4*YU3>>) 
YP«»IMEU4)  »  lYlie)  ♦  Y  ( 13)  ♦  Y  ( 15)  +  Y  (  17)  -  <4»Y<14>>) 

YPRIME  ( 13  >  «  (YIU)  ♦  Y  ( 14  )  ♦  Y  ( 16)  ♦  Y  ( 18 )  -  <4*YU5))> 

YPRIME  U6)  »  ( Y  <  12)  +  I2IYIIS))  ♦  YU9)  -  (4*Y(16>>> 

YPRIME  (17)  -  <(2*Y(14))  ♦  (2*YU8>)  -  (4*Y<17>>> 

YPRIME  (18)  ®  (YUS)  ♦  Y  ( 17)  ♦  Y  (  19 )  ♦  Y(20>  -  (4*YU8))) 

YPRIME  (19)  «  (Y  ( 16)  ♦  (2*Y(18)>  ♦  Y(21>  -  (4*Y(19U> 

YPRIME  (2.0)  »  ((2*Y(10I)  ♦  (2IYI21))  -  <4*Y(20)>> 

YPRIME (21 )  =  Y ( 19 )  ♦  (2*Y(20))  -  C*Y(21)I 

RETURN 

END 
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SUBROUTINE  FCN1  IN, X , Y , YPRIME ) 
INTEGER  N 

REAL  Y  N )  ,  YFR I  ME  !  N )  ,  X 


YPRIME  < I )  =  0.0 
YPRIME ( 2 )  =  Y  <  3 )  ♦  2.*Y!1>  ♦  YI12)  -  4.*Y(2) 
YPRIME  <  3)  =Y  (  4MY!  2)  *Y  1  i4>  -4.  *Y  (  3UYI1I 
YPRIME !  4 ) “Y  <  5MY!  3) +Y < 13) -4 . #Y (  4MY<1> 
YPRIME (  3) =Y  <  6 ) ♦Y !  4 ) ♦ Y < 16 ) -4 . *Y <  3)*Y(I> 
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♦  2.#Y(18> 

YPRIME (24) 

a 

C3K  Y ( 25 ) 

♦ 

C3B*Y(23>  -  4 . *Y (24 ) 

♦  2 .  itY  (  19 ) 

YPRIME (23) 

m 

C4*Y (26) 

♦ 

C4B*Y(24>  -  4.*Y(25> 

♦  2.*Y(20) 

YPRIME (26) 

a 

2.*Y(23> 

- 

4.*Y<26)  *  2.*Y(21> 

^4] 

r « 

0 


S:'.n 
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